function [cs, cv, cl, cc] = splitclustercoords(c, v, k, d)
% splitclustercoords  - split coords of one cluster to subclusters
%
% FORMAT:       [cs, cv, cl, cc] = splitclustercoords(c, v [, k [, d]])
%
% Input fields:
%
%       c           coordinates of values
%       v           values
%       k           size threshold for sub-clusters (default: 3)
%       d           optional 1x3 minimum distance [default: [2, 2, 2]]
%
% Output fields:
%
%       cs          list of cluster sizes
%       cv          Cx1 cell array with lists of values
%       cl          Vx4 lister of cluster voxels
%       cc          Cx1 cell array with lists of coordinates

% Version:  v0.7f
% Build:    8110521
% Date:     Nov-05 2008, 9:00 PM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% check arguments
if nargin < 2 || ...
   ~isa(c, 'double') || ...
    ndims(c) ~= 2 || ...
    size(c, 2) ~= 3 || ...
    size(c, 1) < 4 || ...
    any(isinf(c(:)) | isnan(c(:)) | c(:) < 1 | c(:) > 512 | c(:) ~= round(c(:))) || ...
   ~isa(v, 'double') || ...
    numel(v) ~= size(c, 1) || ...
    numel(v) ~= max(size(v)) || ...
    any(isinf(v) | isnan(v))
    error( ...
        'BVQXtools:BadArgument', ...
        'Bad or missing argument.' ...
    );
end
nc = size(c, 1);

% linearize given values
v = v(:);

% make sure coordinates are as low as possible (for clustering volume)
mnc = min(c) - 1;
if any(mnc > 0)
    c = c - mnc(ones(1, nc), :);
end
mxc = max(c);

% get minimum value and patch values to be >= 1
mnv = min(v);
mv = (1 - mnv) + v;

% create volume filled with a value that is definitely lower than anything
clvn = -(abs(mnv) * 2 + nc);
clv = zeros(mxc) + clvn;

% fill with patched values
clv(sub2ind(mxc, c(:, 1), c(:, 2), c(:, 3))) = mv;

% check there is only one cluster
numclusters = clustercoordsc(clv > 0, 3, 1);
if numel(numclusters) > 1
    error( ...
        'BVQXtools:InternalError', ...
        'Subclusters only valid for single cluster.' ...
    );
end

% what threshold for subclusters
if nargin < 3 || ...
   ~isa(k, 'double') || ...
    numel(k) ~= 1 || ...
    isinf(k) || ...
    isnan(k) || ...
    k < 1
    k = 3;
elseif k >= (nc ^ (2 / 3))
    k = round(nc ^ (2 / 3));
else
    k = round(k);
end

% what distance %% TO BE IMPLEMENTED %%
if nargin < 4 || ...
   ~isa(d, 'double') || ...
    numel(d) ~= 3 || ...
    any(isinf(d(:)') | isnan(d(:)') | d(:)' < 1.5 | d(:)' > (mxc / 2))
    d = [2, 2, 2];
end
dmid = 1 + ceil(d);
dsiz = dmid + ceil(d);
dmsk = false(dsiz);
for xc = 1:dsiz(1)
    for yc = 1:dsiz(2)
        for zc =1:dsiz(3)
            dd = ([xc, yc, zc] - dmid) ./ d;
            if sum(dd .* dd) < 1
                dmsk(xc, yc, zc) = true;
            end
        end
    end
end
[dmsi{1:3}] = ind2sub(dsiz, find(dmsk(:)));
dmsi{1} = dmsi{1} - dmid(1);
dmsi{2} = dmsi{2} - dmid(2);
dmsi{3} = dmsi{3} - dmid(3);

% prepare other arrays
ci = zeros(1, nc);

% sort values and coordinates
[mv, mi] = sort(mv, 'descend');
v = v(mi);
c = c(mi, :);

% set max peak + distance to -1
ci(1) = 1;
cs = zeros(1, nc);
tcrd = validcoords( ...
    [dmsi{1} + c(1, 1), dmsi{2} + c(1, 2), dmsi{3} + c(1, 3)], mxc);
clv(tcrd) = -1;

% iterate over remaining values and
numclus = 1;
for rc = 2:nc
    
    % get coordinate
    dc = c(rc, :);
    
    % check for other clusters
    df = max(dc - 1, 1);
    dt = min(dc + 1, mxc);
    cltouch = lsqueeze(clv(df(1):dt(1), df(2):dt(2), df(3):dt(3)));
    cltouch = max(cltouch(cltouch < 0));
    
    % touches cluster -> mark to one with highest peak
    if cltouch > clvn;
        clid = -cltouch;
        ci(rc) = clid;
        clv(dc(1), dc(2), dc(3)) = -clid;
        cs(clid) = cs(clid) + 1;
        
    % does not touch
    else
        
        % increase number of identified sub-clusters
        numclus = numclus + 1;
        cs(numclus) = 1;
        tcrd = validcoords( ...
            [dmsi{1} + dc(1), dmsi{2} + dc(2), dmsi{3} + dc(3)], mxc);
        clv(tcrd(clv(tcrd) >= 0)) = -numclus;
        ci(rc) = numclus;
    end
end
cs = cs(1:numclus);

% apply sub-clustering k-threshold
for rc = numclus:-1:1
    
    % smaller ?
    if cs(rc) < k
        
        % for each voxel find adjacent cluster with highest number
        cci = find(ci == rc);
        ccoords = c(cci, :);
        for scc = cs(rc):-1:1
            % get coordinate
            dc = ccoords(scc, :);

            % find other clusters
            df = max(dc - 1, 1);
            dt = min(dc + 1, mxc);
            cltouch = lsqueeze(clv(df(1):dt(1), df(2):dt(2), df(3):dt(3)));
            cltouch(cltouch == -rc | cltouch < -nc | cltouch >= 0) = [];
            if isempty(cltouch)
                clid = 0;
            else
                clid = -max(cltouch(cltouch < 0));
                if clid <= numel(cs)
                    cs(clid) = cs(clid) + 1;
                else
                    clid = 0;
                end
            end
            
            % reasign
            clv(dc(1), dc(2), dc(3)) = -clid;
            ci(cci(scc)) = clid;
        end
        
        % then remove cluster
        cs(rc) = [];
        numclus = numclus - 1;
    end
end

% renumber clusters
cn = find(ci == 0);
c(cn, :) = [];
ci(cn) = [];
v(cn) = [];
clnum = unique(ci);
clnid = zeros(1, max(clnum));
clnid(clnum) = 1:numel(clnum);
ci = clnid(ci);

% add coordinate
if any(mnc > 0)
    c = c + mnc(ones(1, size(c, 1)), :);
end

% create output
cv = cell(numclus, 1);
cl = [c, ci(:)];
cl(:, end) = ci(:);
cc = cell(numclus, 1);
for rc = 1:numclus
    cv{rc} = lsqueeze(v(ci == rc));
    cc{rc} = c(ci == rc, :);
end


%%% sub functions

% validcoords
function c = validcoords(c, m)
c(any(c < 1, 2) | any(c > m(ones(1, size(c, 1)), :), 2), :) = [];
c = sub2ind(m, c(:, 1), c(:, 2), c(:, 3));
